This R notebook conducts various tests using the method outlined in Toward automated extraction and characterization of scaling regions in dynamical systems.
# Dependencies
library(ggplot2)
library(LaplacesDemon)
source("./generic.R")
source("./deshMethod.R")
states_visited <- c()
acfs <- c()
# Specific to autoregression
a = 0.95
b = 2
exp_m <- 1 + 2*(a/(1-a))
steps = 10000
get_next_step <- function(prev) {
w <- rnorm(n=1, mean=0, sd=1)
next_val <- a * prev + b * w
return(next_val)
}
get_init_state <- function( ) {
init_sd <- b^2 / (1-a^2)
init_state <- rnorm(n=1, mean=0, sd=init_sd)
return(init_state)
}
# example method with one simulation using one set of p, q
# also plots histograms here.
n <- 10 # chosen
p <- 3 # chosen
q <- 3 # chosen
steps = 10000
a <- 0.9
r <- 0
states_visited <- asim(get_init_state(), get_next_step, steps)
acf_df <- c_est_series(states_visited)
# limiting the search to a finite number, since takes a long time computationally
stop_search <- which(acf_df$acf <= 0)[1]
stop_search <- 1.5 * stop_search
y <- get_log(acf_df$acf)[1:stop_search]
x <- c(0:(stop_search-1))
lim <- length(y)
results_ar1 <- desh_method(x,y, lim, n)
[1] 58
Time difference of 0.2918248 secs
new_weights <- get_desh_weights(results_ar1, p, q, r)
[1] 1128
vars_ar1 <- get_slopes_sides(results_ar1, new_weights, TRUE)
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
vars_ar1
exp_m <- 1 + 2*(a/(1-a))
sprintf("Resulting RHS : %f, expected time: %f", vars_ar1$best_rhs, exp_m)
[1] "Resulting RHS : 20.932509, expected time: 19.000000"
This is the big test– does the method for many different a’s, p’s, and q’s.
# Generates simulations for each possible a we want to test, storing results in a dictionary.
#poss_as <- c(-0.9, -0.7, -0.5, 0.5, 0.7, 0.9, 0.99)
poss_as <- c(0.5, 0.7, 0.9, 0.99)
steps = 10000
dict <- c()
for(a in poss_as) {
states_visited <- asim(get_init_state(), get_next_step, steps)
acf_df <- c_est_series(states_visited)
dict[as.character(a)] <- acf_df
}
# Then we try the method on each simulation.
# Start by finding the autocorrelation functions.
poss_ps <- c(1, 2, 3)
poss_qs <- c(1, 2, 3)
poss_rs <- c(0)
old_ms <- c()
new_ms <- c()
ps <- c()
qs <- c()
as <- c()
real_m <- c()
results_ar1s <- list()
for(a in poss_as) {
acf <- dict[as.character(a)]
acf_df <- data.frame(acf)
colnames(acf_df) <- c("acf")
stop_search <- which(acf_df$acf <= 0)[1]
stop_search <- 2 * stop_search
if(a < 0) {
n <- 3
} else if(stop_search > 500) {
stop_search <- 500
} else {
n <- 10
}
if(stop_search < 5) {
stop_search <- 10
}
y <- get_log(acf_df$acf)[1:stop_search]
x <- c(0:(stop_search-1))
lim <- length(y)
results_ar1 <- desh_method(x,y, lim, n)
results_ar1s <- append(results_ar1s, results_ar1)
}
[1] 20
Time difference of 0.0153501 secs
[1] 36
Time difference of 0.1190212 secs
[1] 160
Time difference of 6.823408 secs
[1] 500
Time difference of 8.277671 mins
# Then try the method on the given autocorrelation functions.
old_ms <- c()
new_ms <- c()
ps <- c()
qs <- c()
as <- c()
real_m <- c()
rs <- c()
for(a in poss_as) {
print(a)
start_time <- Sys.time()
c <- 4
M_upper_bound <- 6000
triangle_df <- find_m(acf_df, c, M_upper_bound)
m_suggest <- triangle_df[triangle_df$poss_Ms >= triangle_df$ct_ints,]
m_suggest <- m_suggest$poss_Ms[1]
m_suggest <- t_int(m_suggest, acf_df)
acf <- dict[as.character(a)]
acf_df <- data.frame(acf)
colnames(acf_df) <- c("acf")
stop_search <- which(acf_df$acf <= 0)[1]
stop_search <- 1.5 * stop_search
if(a < 0) {
n <- 3
} else if(stop_search > 500) {
stop_search <- 500
} #else {
#n <- 10
#}
if(stop_search < 5 ) {
stop_search <- 10
}
y <- get_log(acf_df$acf)[1:stop_search]
x <- c(0:(stop_search-1))
lim <- length(y)
results_ar1 <- desh_method(x,y, lim, n)
for(p in poss_ps) {
for(q in poss_qs) {
for(r in poss_rs) {
exp_m <- 1 + 2*(a/(1-a))
real_m <- append(real_m, exp_m)
old_ms <- append(old_ms, m_suggest)
new_weights <- get_desh_weights(results_ar1, p, q, r)
vars_ar1 <- get_slopes_sides(results_ar1, new_weights, plot_title=sprintf("a:%f, p:%f, q:%f", a, p, q)) #automatically plots RHS
new_ms <- append(new_ms, vars_ar1$best_rhs)
ps <- append(ps, p)
qs <- append(qs, q)
as <- append(as, a)
rs <- append(rs, r)
}
}
}
end_time <- Sys.time()
print(end_time - start_time)
}
[1] 0.5
[1] 15
Time difference of 0.003687143 secs
[1] 10
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 10
[1] 10
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 10
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 10
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 10
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 10
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 10
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 10
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
Time difference of 59.63924 secs
[1] 0.7
[1] 27
Time difference of 0.04512501 secs
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 136
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
Time difference of 59.47326 secs
[1] 0.9
[1] 120
Time difference of 2.51026 secs
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 5995
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
Time difference of 1.037575 mins
[1] 0.99
[1] 500
Time difference of 7.548814 mins
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
[1] 119805
Warning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true densityWarning: sum(weights) != 1 -- will not get true density
Time difference of 11.9294 mins
# Formatting and storing results
results <- data.frame(old_ms, new_ms, ps, qs, as, real_m)
colnames(results) <- c("Old Method", "New method", "p", "q", "a", "Real Time")
write.csv(results, "./results.csv", row.names = FALSE)
results
# plotting the densities for the best p, q for each a
for(a in poss_as) {
results_a <- results[results$a == a, ]
results_a$err_vs_real <- (results_a["Real Time"] - results_a["New method"]) / results_a["Real Time"]
selected_vals <- results_a[which.min(abs(results_a$err_vs_real$`Real Time`)),]
acf <- dict[as.character(a)]
acf_df <- data.frame(acf)
colnames(acf_df) <- c("acf")
stop_search <- which(acf_df$acf <= 0)[1]
stop_search <- 1.5 * stop_search
if(a < 0) {
n <- 3
} else if(stop_search > 500) {
stop_search <- 500
} else {
n <- 10
}
y <- get_log(acf_df$acf)[1:stop_search]
x <- c(0:(stop_search-1))
lim <- length(y)
results_ar1 <- desh_method(x,y, lim, n)
print("yo")
print(results_ar1)
new_weights <- get_desh_weights(results_ar1, selected_vals$p, selected_vals$q)
vars_ar1 <- get_slopes_sides(results_ar1, new_weights, plot=TRUE)
}
[1] 3
Time difference of 0.0002748966 secs
[1] "yo"
[1] 0
Error in density.default(results$slopes, bw = "nrd", adjust = 0.25, kernel = "gaussian", :
argument 'x' must be numeric